Eleonora De Leonardis, Benjamin Lutz, Sebastian Ratz, Simona Cocco, Remi Monasson, Alexander Schug, and Martin Weigt
尽管非编码RNA具有重要的生物学意义,但它们的结构描述仍然具有挑战性。通过使用急速增长的序列数据库,我们用直接耦合分析(Direct-Coupling Analysis,DCA)的方法分析同源序列上的核苷酸协同进化,来检测核苷酸之间的关联。在展示的核糖开关数据集中,我们发现,DCA结合一般的Nussinov算法能系统的提高RNA二级结构的预测结果,对比于传统的基于相互作用信息(Mutual Information, MI)协方差方法。更重要的是,我们发现,在三级结构关联上DCA的结果也增强。通过整合这些预测结果到分子模型工具,对比仅仅只是用二级结构信息,能系统的提高三级结构的预测结果。
实验工作和基因组序列分析揭示了在细胞中RNA有广泛的作用[1-5]。除了遗传信息的传输,非编码RNA能催化生物化学反应,并在多种调节过程中有关键作用。还有一些功能RNA在它们的单链信息或是RNA-蛋白质复合物环境下表现很重要,另一方面,功能是与RNA的三维结构直接相关[6]。获取这些结构知识对理解其功能非常重要。然而,通过实验测定RNA的结构仍然具有挑战性。在RCSB的蛋白质数据库(Protein Data Bank,PDB)[7]里的全部结构数据中,只有少于6%的数据包含RNA。由于测序技术的进步,许多在不同生物体里RNA已经被测序,并分类到同源家族的Rfam数据库中[8]。然而,在目前列出的2450个家族里,只有59(2.4%)个包含至少一个PDB结构数据,然而却有954(39%)个包括超过100条序列,甚至566(23%)个超过200条序列。因此为了充分使用海量的序列信息,计算机辅助结构预测方法是一个很有前景的互补技术。不幸的是,目前最先进的计算机辅助方法在未知RNA结构预测上也很难使RMSD(均方根差)值降到8-12埃[9-10],远高于能达到2-3埃精度的X射线晶体结构。通常情况下,计算机辅助方法[11-13]局限于众多限制,包括结果不定,或是要求专家参与,复杂的拓扑结构不能够可靠预测[14-15]。一些在此方向上获得的最好的结果是通过结合计算机辅助和实验方法:突变研究得到的跨残基关联作为约束加入到计算模型[16-17]。
在相关的蛋白质结构预测领域,近年通过探究残基协同进化,在残基关联预测方面有了显著的进步[18]。这些新方法基于两种思路:(i) 两残基(甚至可以在一级序列上相距很远)间的三级结构关联导致氨基酸占位的相关模式,这个可以被多序列比对(Multiple Sequence Alignments,MSA)统计分析方法检测到。(ii) 本地相关性测量,比如相互作用信息(MI),因为传递效应而令人困惑:即使不直接连接而是通过相邻联系到一个共同的中间残基的两个位置也将相关。像DCA[19,20]和与之类似的方法[21,22]它们的优点是能够区分直接和间接效果来推断间接相关的直接耦合作用。这导致在预测关联数据的精度大幅提高,可用于预测三级和四级蛋白质结构,例如球蛋白[23,24],蛋白质复合物[25-27],积极构象[28]或膜蛋白[29]。
我们的目标是提出一个有效的基于协同进化分析现有的同源序列的RNA二级结构和三级结构预测流程。用协方差模型进行比较RNA序列分析是一致公认的[30,31]:MI能够成功的用来推断基本配对和预测二级结构[30,32,33]。尽管已经观察的一些三级结构间的联系表明MI的重要性[32,33],然而,有人认为,在三级结构中存在关联的非规范碱基对较少表现出协变性[35]。这里我们发现,在将本地协方差分析方法换为DCA时,这些微弱的信号却能在三维结构关联上显著增强。使用经过严格筛选的具有复杂结构的核糖开关家族数据集,我们发现DCA能够有效的整合进现有的RNA二级结构和三级结构预测的工具中。我们表明,结合类似于Nussinov算法之类的标准方法[36],相比基于MI的预测方法DCA能系统的提高二级结构的预测程度。将DCA整合到Rosetta中[37],一个非常成功的RNA结构预测工具[10],我们提出了一个完全自动化的三级结构预测方案。相对于仅仅使用Rosetta,其结果有显著提高,并且比那些RNA-puzzle实验中包含的工作流更具有竞争力,因为那些方法部分的需要额外的实验信息或者是专家的操作[9,10]。
当运用到同源家族的RNA上时,DCA能够检测二级结构和三级结构关联,这反过来又可以被集成到现有的用于RNA结构预测工具。为了这样做,我们已经开发出一种专门的流程,如图1所示。从简单的RNA序列开始:
此操作流程是完全的模块化,可以在每个模块中改进——更好的序列比对,更精确的DCA信息提取,整合其他(例如实验方面)的知识到分子建模等等。这些都可以方便的实现,并会改善总体性能。
为了测试基本的操作流程,我们确定了测试家族数据集。特别的我们集中在核糖开关上,它具有一定的复杂结构。为了产生一个代表性的RNA家族数据集在我们的研究中用于基因组分析和结构预测,我们从Rfam数据库中通过以下标准系统地选择了核糖开关家族:
如果有多个满足这些条件的PDB结构,我们选择最高精确度的那个。我们发现有六个家族(见附录表S1)满足这些条件,那些不满足条件的家族则不包含在我们的分析中。有趣的是,这些核糖开关是很多研究常用的目标[38,41]。
对于这些家族,我们从Rfam数据库中提取了多序列比对和一致的二级结构[42]。为了适用于特定的RNA序列的结构预测,通过去除关联精选那些一致的二级结构,这在特定序列会不满足碱基基本配对,并有可能通过添加邻近二级结构的基本配对延伸螺旋链。见附录表S2。
这些选择标准不应被理解为我们的方法的适用性限制。正如附录图SI1所示,尽管使用了至少250个子样品,DCA提取的关联信息也不会改变多少。整个操作流程同样适用于超过100个核苷酸的序列,只是在分子建模步骤的计算成本有所增加。
为了证明我们的操作流程除了测试集外仍具有广泛的适用性,我们还分析了两个未知3D结构的RNA家族,即glnA核糖开关家族RF01739和C4反义RNA家族RF01695。我们提供了未知结构预测的十个最高的打分聚类,见附录数据(pdb文件),可以进行实验上的检测。
接下来我们简单的回顾DCA的主要部分(更多详细的描述包括具体技术细节请见附录,其中主要步骤[19]是总结和适用于RNA的特异性):DCA旨在通过广义Potts模型(或者说马尔可夫随机场),MSA作为输入,来模拟序列的变异性。在L个核苷酸或间隔中,它分配一个概率给每个比对序列,
其中,指的是位于
处的核苷酸
和位于
处核苷酸
之间的直接耦合,并且
是局部偏差(场),仅涉及处于位置
的单个核苷酸。分母Z起到归一化作用。参数必须调整以便再现从MSA获取的核苷酸实验统计数据:
对于所有位置和
以及所有的核苷酸
和
。这里,
贡献核苷酸
在位置
出现的频率,
指的是序列在位置
核苷酸为
并且在位置
核苷酸为
的概率。使用平均场方法[19]来估计参数。为了能够根据它们的耦合能力给配对位置排序,5x5维度的矩阵
压缩到一个量化的耦合打分
,详情见附录。对于蛋白质,研究表明此过程十分优于像MI这种局部相关测量,
其子序列平均值与修正值相乘(即MIapc)[43]。
RNA分子的二级结构的协同进化已经众所周知[44];沃森-克里克碱基互补配对需要配对核苷酸在进化过程中协调突变。它是目前为止部分上最先进的预测RNA二级结构[45]和RNA多序列比对方法[34]。大多数预测二级结构的方法受到Zuker的最小自由能算法的启发[46]并且可以整合来自残基共变的附加的信息[47];这些方法的可靠性关键取决于实验确定的热力学参数。因此我们只关心评估使用的统计序列信息,我们遵循一套简单的流程。其基于Nussinov算法[36],原本设计用来通过动态规划编程发现二级结构碱基最大互补配对数,其一般形式由Eddy和Durbin[30]提出,包含了MI作为残基共变测量。
我们也使用了一般Nussinov算法,但是其中的MI替换为其它的多种打分,包括MI本身,也还有MIapc和DCA打分Fapc[19,20]。最小的发卡环长度设置为3个核苷酸。注意这些算法并不作为二级结构预测,但是它们的目标是评估在其他方面相同的算法,但使用DCA打分而不是基于MI打分的影响。将DCA打分整合到类似RNAalifold[47]这种方法中将会非常有意思。当然,这不在本片文章讨论范围内。
Nussinov算法需要一个关联打分矩阵作为输入,Durbin提出的一个变体仅仅使用MI,由于其值均为正数,会导致过于配对。除此之外,映射到特定的RNA序列上时,可能其配对的核苷酸表示出强烈的协变性,但是形成的不是互补的沃森-克里克配对或是次稳定配对。为了避免这两种情况,我们对于特定的目标序列采取下面的策略来构造配对打分矩阵:
这些规则保证了矩阵上最多有L个正的入口,因此避免了过配对。得到的预测的二级结构与之对应的PDB上的二级结构进行对比(评估所有来自PDB文件和Rfam比对的残基,即,我们可以得到预测结果并进行评估)。
二级结构预测的精确度与保留的协变打分数目之间的关系在附录中进行了讨论:一个潜在的问题是实际的碱基对的数目比例是序列的长度,还是
。对于用DCA预测,Rfam一致结构的敏感度和精确度(使用PDB二级结构评估)在截断L处趋于相同,此外在相同敏感度的情况下使用MI和MIapc预测的结果精确度相对较低,详情见附录图SI2。
为了在我们的预测中表现三级结构关联增强,我们采取了下面的流程:
在全部的前个预测结果中,局部
值比全局
值是一个更加严格的增强标准,因为TP比值凭经验能观察到几乎单调降低(见结果)。接下来,可视窗口大小
选为全部元素的10%,此值是
值在局部的分辨率(小
)和可靠性(大
)之间的折衷。
对于三级结构预测,我们遵循一般Rosetta的过程[48]。首先,基于二级结构信息搭建理想化的双螺旋结构。第二步,RNA连接点和环结构基元加到二级结构上。这些元素组合成完整的模型,同时考虑预测的三级关联作为能量的约束。对于这些约束,基于[49]将两残基间预测的残基-残基关联映射到原子间的距离约束。基于这些约束和其内能模型,Rosetta生成上千种结构。这些预测的结构按照它们的打分排序,然后以4埃为临界值进行聚类分析来确定具有代表性的结构。详细的过程描述见附录。
首先在RNA家族上测试DCA的表现是用直接耦合打分作为Nussinov动态规划算法的二级结构预测的输入,详情见材料与方法。可能有人会说Watson-Crick碱基配对的核苷酸之间有非常强且直接耦合的关系,所以二级结构预测不一定需要进行DCA。如图2所示,我们发现,尽管如此,在相同的精确度(即真阳性预结果在所有预测结果的比例)的情况下敏感度(即有更多的二级结构关联被预测出来)有很明显提高(真阳性(TP),假阳性(FP)和假阴性(FN)等相关术语的定义见图题),与之对比的是使用MI。你也可以在蛋白质结构预测中发现,通过应用平均积修正改进[43]可以提高使用MI预测结构的结果,但是这种提高远远小于通过DCA方法获取的。在图3中,你可以看到分别使用DCA和MI预测二级结构的结果以及基于PDB文件中三维结构信息的参考结构(见附录图SI3关于MIapc的对比)。需要注意的是两者都有的错误,两者的假阳性(由绿线连接的配对)和假阴性(蓝色填充的基本配对)大部分出现在靠近晶体结构序列和Rfam比对序列的非比对区域(灰色颜色标记)间。对于这些区域,我们无法正确定义Nussinov打分矩阵,因此错误更可能发生。
此外,相对于给定的实验结构,一些明显的DCA错误预测可能会指向特定的PDB序列比对错误。一个特别的例子是2gdi的第三个支干(基本配对15-25,16-24,17-23,以及18到22的环)。首先,假阳性位于一个高度间隔区域:在比对中只有59%的序列不包含任何在这3个基本配对所占的6个位置上的间隔。另外分析表明,在无间隔序列中,大多数(54%)是完全与DCA预测的基本配对相同的,只有2%是与PDB结构完全相同。31%的序列与DCA结果和PDB结构都相同,其他所有的序列在二级结构上部分相同。在这种情况下,两个核苷酸如果能够形成Watson-Crick配对或是wobble配对就称为碱基配对相同。
当我们分析RNA家族的大序列比对时我们能够提取二级结构之下的关联信息吗?有人认为,三级结构的核苷酸之间关联的协同进化信号比二级结构信号弱的多[35],并且在许多RNA家族中可能与噪音无法区分。
我们注意到这幅图是将局部协同进化方法(MIapc)替换为全局方法(DCA)。图4中的曲线展现了经过了六个RNA家族的平均,真阳性预测结果与总的阳性结果的数量无关,见附录图SI4有每个家族的详细结果。上面的两幅图使用了严格的大小为4埃的关联截断,而下面的图则使用相对宽松的8埃最为临界值,相比所选的代表性结构,其目的是反映出RNA家族成员里RNA结构的可变性。
第一眼看上去,使用MI,MIapc和DCA得到的结果看上去非常相近,见左边的两幅图。这导致一个结果,即最强的信号实际上是由二级结构关联决定,见底下两幅图中的亮线,这条线表明三级结构关联中得分最高的预测,其在一开始贴近于0。为了能更仔细的评估DCA在预测三级结构关联的表现,我们因此去掉了二级结构和相邻配对(对于每一个二级结构配对我们同样从所有的预测数据集中去除掉
)。正如图4右边所示,相比之于MI和MIapc,剩下的DCA预测结果包含基本上更高的关联部分,但是他们仍然还是离黑线代表的最好的预测结果较远,与此同时列出所有非关联之前的关联信息。
当处理大量的RNA家族数据集时,一个公认的DCA相比MI和MIapc的改进也被发现。附录图SI5展示了15个家族的结果,这些家族有超过1000条序列,精度好于3埃的PDB结构以及较丰富的三级关联(包括简单发卡环),但是没有对序列长度或功能类的限制。大型核糖体RNA也包括在统计数据中,因为它们在序列长度上和数量上与其它所有RNA家族都不同。
Rosetta是基于二级结构信息(SSI)来搭建理想的双螺旋模型,然后基于片段数据库增加环形区域到建议的上千个结构单元来组成一个大的序列部分。所有的片段通过蒙特卡罗过程使得内能最小然后组合在一起。核糖核酸之间用DCA或者是MI预测出的三级残基间关联可以作为基于原子的约束[49]来修改能量部分,即它们通过偏置打分能量引导结构预测。接下来,我们关心的是六组不同的结构约束/信息集合:(I)在RFAM的MSA中的基于SSI的一致二级结构,(II/III)SSI加上25个或具有代表性的100个最高的MIapc排名的核苷酸配对,(IV/V)SSI加上25个或者具有代表性的100个最高的DCA排名的核苷酸配对,(VI)SSI加上全部的三级内环残基关联映射。集合I和VI定义了最坏或最好的情况下的参考值(即没有对全部的三级关联信息)用于基于Rosetta能量打分的RNA结构预测[50]。他们表现出了改进的可能性,可以通过增加三级结构约束到SSI中来实现。集合II到V提供了,用于预测关联信息的MIapc和DCA,这两种可供选择的策略使用了很多但低质量的关联预测,与之对比的是使用相对较少但高质量的关联预测。虽然第一组可能会由于假关联预测导致结构预测为错误结构,第二组可能会遗漏掉那些本身低协同进化信号的三级关联的聚类,参见附录图SI6中的关联映射。
我们通过Rosetta将这些内部残基关联预测作为约束增加到RNA结构预测中去[37],遵循[48]里面的实验报告。
我们按照这个实验报告对所有的六个核糖开关例子进行处理。表1列出了那些最好的Rosetta分数的结构模型的RMSD值,以及在10个Rosetta分数里具有代表性的5个最小的RMSD,假定附加信息(例如已知的实验约束[51])将允许选择最准确的预测。除了一个案例,在仅仅提供SSI(集合I)的情况下增加三级关联(集合VI)强烈的改善了RMSD值,显示了Rosetta用于进行分子建模的外部三级关联信息的值。
一般来看DCA指导的预测结果优于SSI和MIapc指导的。在某些情况下使用全部关联信息预测的准确度几乎达到此地步,该预测也与天然结构的预测相覆盖,见图5。对于集合II和III,基于作为局部协变测量的MIapc,我们观察到不同的核糖开关的质量参差不齐的预测(在前面最好的10个预测结果,包括集合II和III:平均值 = 11.8埃,最大值 = 16.3埃,最小值 = 7.5埃)。相反,对于集合IV和V,基于DCA,我们观察到对比集合I预测质量有显著提高,在大多数情况下也超过了集合II和III(在前面最好的10个预测结果,包括集合IV和V:平均值 = 9.6埃,最大值 = 12.4埃,最小值 = 6.7埃)。
还有一些特殊的情况值得注意:查看附录图SI6里2gdi的关联映射,MI和DCA预测的结果似乎相同,但是表1所示的预测结构表明使用DCA具有更高的准确度。仔细查看关联映射后原因就变得清晰了。MI预测完全忽略了三个天然态关联的聚类,然而DCA预测则找到了所有的聚类。这表明很少预测可能对最终的预测精度有重大影响,当他们加入正确的距离条件约束到某一区域,而在其他方面不加约束,从而减少了大量可行的三维结构的样本空间。这个发现在预测3owi的相对低的准确度上的到确证,除去一些假阳性。预测的关联映射主要由二级结构决定,还有一些相对有限的三级结构信息也加进我们的方法中。
另一方面,2gis表现出相对大量的假阳性,它与任何天然态的二级结构或三级关联相距甚远并且完全由Rosetta预测出结构。实际上,观察到的RMSD值在我们的测试组里至少稳定。
另一个有趣的例子是3vrs,它暴露出了Rosetta建模中的一个问题。其结果显示DCA优于全关联映射。我们发现全关联映射预测实际上有更好的结果(约9埃),但是它们不被打分系统检测到,你可以见附录图SI7。聚类的最低分数始终是12到16埃RMSD。更近一步对Rosetta准确度的限制甚至在全关联映射结果的情况下,事实上我们仅仅提供天然态的残基-残基关联列表,但没有关于天然态原子间距离的信息。
更多信息可以见附录资料。
这些测试表明,除开弱协同进化的三级结构引起的信号联系,DCA结果能够很好地集成到RNA 3D结构预测系统,并且有效地提高预测精度。对于一个给定的校准(Stockholm或Fasta格式)在标准的台式电脑下获取DCA预测最多只需要几分钟。DCA流程基于三个步骤,参见补充数据:在L是序列长度,M是序列编码的条件下,权重是的内存空间和
的时间复杂度,估算协方差是
的内存空间和时间,矩阵求逆是
的内存和
的时间。映射残基-残基的接触到原子-原子的接触只需要几秒钟(
时间和空间需求之下)。Rosetta则需要最大的计算机工作量。第一步是生成循环区域模型。对于每一个检查过的核糖开关,每个循环区域约4000个模型的生成需要大约3CPU天。对于Rosetta的最后一步,从螺旋循环区域部分装配完整的核糖开关,我们花了36个CPU天让每个核糖开关达到了足够的采样(2000-5000个模型)。整个流程完全自动化,并且在任何一个步骤中不依赖任何人类接触预测或装配。
这种方法直接适用于本文提到的之外的其他RNA的分析。为了展示这一点,我们提供了RNA家族的两个盲测:第一个RNA家族是谷氨酰胺合成酶基因主题(RF01739),即谷氨酸盐粘合物核糖开关。第二个例子我们取缔了核糖开关。C4反义RNA家族(RF01659)中发现了噬菌体感染细菌。这两个例子都符合我们之前选择标准的序号和长度,但没有实验性PDB结构是一致的。每个RNA家族的特定的成员序列必须选择Rosetta建模:为了将定位错误的风险降到最低,我们选择序列特定的DCA预测和Rfam一致二级结构。作为结果的PDB模型(每个家族十个最佳评判模型,使用SSI和100DCA预测构筑)都包含在补充数据中,2D和3D的表达如图6所示,接触图示的补充如图SI8所示。
最后,我们想提出一些关于如何提高预测的未来发展方向。未来工作中肯定会出现预测三级结构的改良版Rosetta协议和(或)评分方案,而我们的重心放在了使用DCA统计预测时出现的合理改进。
首先,通过多重序列比对确定RNA是一个复杂且未完全解决的问题,同时用于生成大型Rfam的测量分析系统程序可能仍然缺乏准确性。结果中我们已经展示出了一个因比对错误限制了二级结构预测的准确性的例子。为了调查关于三级结构预测中比对质量带来的影响,我们完成了以下数值实验:对于Rfam家族RF00010(细菌核糖核酸酶P类甲级),我们在三种不同测量系统分析下运行DCA:(I)全部Rfam测量系统分析包含6397种序列;(II)340种序列(52)的优质结构对齐序列;(III) sub-MSA Rfam MSA包含序列中包含的结构排列,参见图7。毫无意外的是最小和品质最有限的队列(III)执行结果最差。然而几乎所有在此比对比较优势都通过提高在一个固定数量的序列上的对齐质量或通过提高Rfam完整比对所RNA家族采样对齐而获得的。
即使对于一个既定的MSA,协同进化分析预测也可能在下述综合方法的启发下而得到提升,如[53,54],使用DCA输出(或非常相似的PSICOV[22])时,协同其他的功能,如二级结构预测信息和通过机器学习工具输入溶剂可接触性。此类监督方法若要获得成功,可能受限于RNA在由于少数的Rfam家族具有已知结构典型的情况之下。然而有证据表明,原生RNA结构中接触的实质部分显示了某种程度上的协同进化。为了定量评估这个想法,我们计算出了在协同进化预测中接触富集性的假定值。更准确地说,我们比较了在滑动窗口下原始TP率和随机预测的TP率(适用于不在更高DCA队列中的所有配对,也就是仍能被预测的所有接触,参见材料与方法)。图8表示了当假定值第一次超过意义阈值为0.01的时候,灵敏度预测和精度预测的不同的协方差方法。只有18%的接触包含在内时,MI发生(平均每六个测试家族),而且TP/(TP+FP)只有10%的精度,DCA发现所有接触中的54%只略增加了15%的精度。此外,APC应用于MI的结果在简单MI之上,但结果(灵敏度25%,精度13%)仍然远低于DCA。一个完全随机预测会导致P/(P+N)=5.7%的精度。因此,即使与MIapc相比,DCA在接触中达到了一个更强大的富集度,但是结果精度可能会过低于实际三级结构预测。然而,这值得设置旨在更好地辨识正确预测与错误预测的过滤方法,并进一步提高正确预测的富集度。
表1 基于不同残基关联映射的6个代表性的核糖开关的结构预测
| PDB Rfam | Meff/L | Secondary only | 25 MIapc | 100 MIapc | 25 DCA | 100 DCA | All |
|---|---|---|---|---|---|---|---|
| 1y26 RF00167 | 5.8 | 21.6 19.6 12.6 | 12.3 12.3 12.2 | 12.4 10.2 10.2 | 11.4 8.9 8.9 | 16.2 11.3 11.3 | 7.6 7.6 7.4 |
| 2gdi RF00059 | 32.2 | 16.1 16.1 16.1 | 17.0 16.3 8.3 | 21.5 18.5 17.6 | 7.5 7.5 7.5 | 7.7 6.7 6.7 | 6.3 6.3 6.3 |
| 2gis RF00162 | 10.9 | 20.7 18.4 17.4 | 21.8 21.8 17.5 | 22.2 7.7 7.3 | 22.1 10.5 10.5 | 13.9 13.8 10.6 | 9.9 5.8 5.8 |
| 3irw RF01051 | 11.3 | 17.0 17.0 10.2 | 8.8 7.5 7.5 | 8.5 8.4 7.6 | 9.3 9.2 8.6 | 8.1 7.8 7.8 | 6.5 6.5 6.5 |
| 3owi RF00504 | 20.3 | 13.5 13.5 13.5 | 14.4 14.4 13.3 | 17.1 15.2 14.4 | 14.6 12.4 11.7 | 16.2 12.9 11.6 | 10.0 6.3 6.3 |
| 3vrs RF01734 | 8.3 | 15.8 14.2 12.5 | 15.1 14.6 14.6 | 15.0 14.4 13.4 | 15.6 12.3 12.3 | 11.8 11.8 8.8 | 12.9 12.0 12.0 |
附录文件 可在 NAR Online 上获取。
Agence Nationale de la Recherche project COEVSTAT [ANR-13-BS04-0012-01 to S.C., R.M. and M.W., in part]; Impuls- and Vernetzungsfond of the Helmholtz Association [to B.L., S.R. and A.S.]. B.L. acknowledges hospitality of the UPMC in the initial phase of this work. Funding for open access charge: Deutsche Forschungsgemeinschaft and Open Access Publishing Fund of Karlsruher Institut fur Technologie.